In the case of non-stationary problems, the solution sought depends both on time and space.
In the general strong form, in which the computational problem is described by a partial differential equation
\( \frac{\partial u\left(x,y;t\right)}{\partial t}-{\cal L }\left(u\left(x,y;t\right)\right)=f(x,y;t) \)
the derivative with respect to time now appears
\( \frac{\partial u\left(x,y;t\right)}{\partial t } \)
representing the changes over time of the searched scalar field
\( \Omega\times [0,t] \ni \left(x,t\right)\rightarrow u\left(x,y;t\right)\in R \).
The meaning of the searched field depends on the type of problem being solved.
The type of problem to be solved is "coded" in the operator \( {\cal L }\left(u\left(x,y;t\right)\right) \), describing the modeled physical phenomenon by means of partial differential operators.
Similarly to solving stationary problems (not changing with time), it is necessary to provide boundary conditions (i.e. information about what is happening on the boundary of the simulated area). Additionally, in the case of simulation of non-stationary problems, it is necessary to provide the initial condition (i.e. provide information about the state of the modeled phenomenon at the beginning of the simulation) \( u(x,y;0)=u_0(x,y) \).
There is also a right-hand side in the equation \( f(x,y;t) \) describing the "force" that supplies energy (momentum, mass, etc. depending on the modeled phenomenon) to the system.
There are two methodologies for solving the non-stationary problems. In the first, discretization is performed first with respect to time, and then with respect to space. In the second, discretization is performed first in space, and then over time (this method is called the method of lines). In our considerations, we will use the first method.
Thus, in order to solve the non-stationary problem, we introduce time steps
\( t_0=0 < t_1 < t_2 < \cdots t_{k-1} < t < t_{k+1}< \cdots < t_N \)
and states of the modeled phenomenon in particular time steps, represented by a scalar field \( u(t) \)
\( u_0=u(t_0), u_1=u(t_1), u_2=u(t_2), \cdots, u_{k-1}=u(t_{k-1}), u_k=u(t_k), u_{k+1}=u(t_{k+1}), \cdots, u_N=u(t_N) \)
We approximate the derivative in time using the finite difference method \( \frac{\partial u}{\partial t} \approx \frac{u_{t+1}-u_t}{dt } \)
So we get the equation
\( \frac{u_{t+1}-u_t}{dt} - \mathcal{L}(u) = f \).
Now the question arises, at what moment in time to take the state of the modeled phenomenon, represented by the scalar field
\( u(t) \). We have several possibilities:

1. take the state at the previous time moment \( u_t \),
2. take the state at the actual time moment \( u_{t+1} \),
3. take the state as the combination of the previous and actual time moments \( \alpha u_{t} + (1-\alpha) u_{t+1} \).

The first way gives us the explicit method, also called the Euler method (forward Euler)
\( \frac{u_{t+1}-u_t}{dt} - \mathcal{L}(u_t) = f \).
The second way gives us the implicit method, also called backward Euler
\( \frac{u_{t+1}-u_t}{dt} - \mathcal{L}(u_{t+1 }) = f \).
The third way gives us the implicit method, and depending on how the linear combination is defined, it will be the Crank-Nicolson method, or the alpha method. \( \frac{u_{t+1}-u_t}{dt} - \mathcal{L}(\alpha u_{t} + (1-\alpha) u_{t+1}) = f \).

Explicit methods always allow for a quick solution in linear time, but only small time steps can be used in them (otherwise the simulation will become unstable). Implicit methods allow the use of larger time steps (if they are absolutely stable, and checking if the method is indeed absolutelty stable it requires advanced mathematical knowledge). Implicit methods usually require expensive solvers. In this tutorial we shows how to get linear time solvers for implicit methods on tensor product patches of elements employed by the isogeometric analysis.
These methods are derived from a family of methods called Runge-Kutta methods [1].